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Bethe Ansatz equations and the variational BCS wavefunction. Comparisons 
are provided as a function of coupling strength and electron density. While the 



is sometimes incorrect by an order of magnitude, particularly at half-filling. 
^ \ Finite size effects are also briefly discussed for cases where an exact solution 

I in the thermodynamic limit is not possible. In general, the BCS result for the 

energy gap is poor compared to the exact result. 

1. INTRODUCTION 

The discovery of high-temperature superconductivity has motivated a considerable effort 
over the last decade. In particular, many theoretical models have been proposed to explain 
both the superconductivity and some unusual normal state properties |^. As far as super- 
conductivity is concerned, model building has proceeded on at least two paths. The first 
includes "realistic" model Hamiltonians, often two dimensional (since the CUO2 planes have 
been deemed to be the essential structure for electron conduction). The price of such realism 



is that exact solutions are impossible, and a mean field or BCS 0] solution is required. The 
second path attempts to remedy this deficiency by seeking exact solutions to a restricted 
class of models, often in one (and more recently in infinite) dimension. These models are 
often unrealistic, and are sometimes exactly solvable only in restricted parameter regimes. 
Attempts to bridge the gap between these two modes of theoretical work are becoming more 
commonplace in recent years. Two examples are spin fiuctuation theories, where Fluctua- 
tion Exchange calculations are often compared to Monte Carlo calculations for the Hubbard 
model 1^, and Hubbard- like models that include occupation-dependent hopping rates @]. In 
this latter example BCS calculations have been compared with exact diagonalization studies 
on finite systems. 

In these and other examples, comparison of BCS or some other approximation scheme 
with exact results has always been hampered by some deficiency, often finite size effects (one 
exception is mentioned below). Yet an exact solution to the attractive Hubbard model in 
one dimension has been available for almost thirty years, via Bethe Ansatz techniques 0. 

In this paper a thorough investigation of the attractive Hubbard model will be pre- 
sented, using both Bethe Ansatz techniques and BCS solutions. The study is necessarily 
confined to one dimension. It is not yet clear to what extent conclusions obtained in one 
dimension carry over to higher dimensions. Nonetheless, we were encouraged by the results 
of Quick et al. 0, who performed a study similar to the present one using a model of a 
one-dimensional electron-gas with pairwise-attractive delta-function interactions, first solved 
exactly by Gaudin and Yang [§. This is a one parameter model with a dimensionless 
coupling constant inversely proportional to the electron density. Quick et al. 0] found that 
BCS theory gave an accurate estimate of the ground state energy over all ranges of the cou- 
pling constant. Quantum fiuctuations are most severe in one dimension, so the agreement 
should improve as the dimension is increased. 

Part of the purpose of the present study is to test the accuracy of BCS as "an interpolation 
scheme" [0,|Tni over all coupling strengths. As already noted in Ref. 0, BCS theory is exact 



in the weak and strong coupling limits. The present model has two parameters, electron 



density as well as coupling strength, and the question of how well the interpolation scheme 
works as a function of electron density is addressed. Moreover, the gap to the first excited 
state will also be calculated, as this quantity provides a much more revealing test of the 
BCS approximation than does the ground state energy. 

We begin by outlining the model and the theoretical treatments in the following section. 
For the exact solution we rely on previous work by Lieb and Wu and by Bahder and 
Woynarovich We mention some trivial but important points regarding the solution of 



the resulting equations, for both the ground state energy, and the excitation gap. For com- 
pleteness we also outline the BCS solution, which is straightforward. Results are presented 
in the following section, followed by a summary. 

2. FORMALISM 
A. The Model 

The one dimensional Hubbard Hamiltonian is given by 

H = -tY,{clci^^,a + 4+i,aCia )+ UY,n^1inil, (1) 

where the operator (qo-) creates (annihilates) an electron with spin a at site i and 
is the electron number operator with spin a at site i. The parameter t is the hopping rate 
for electrons and merely sets the overall energy scale. The parameter U is negative in the 
attractive case, so that \U\/t is a dimensionless measure of the coupling strength. The 
second parameter is the electron density, n = Y.i,cr < '<^ia >, where N is the number of 
lattice sites, and the expectation value is taken with respect to the ground state. 

This model has particle-hole symmetry, so we will concern ourselves with electron densi- 
ties less than half-filling. In addition, symmetries exist between the repulsive and attractive 
Hubbard models, such that 

E{M, M'; -\U\) = -M\U\ + E{M, N - M'; \U\) (2) 

where E{M,M'; \U\) is the ground state energy of H (with U = \U\) with M down-spins 



and M' up-spins (M + M' = Ne, the total number of electrons). The total z-component of 
electron spin, S^, is given by = ^{M' — M). For the ground state, Sz = 0. 

B. The Lieb-Wu Equations: Finite Systems 

An exact solution to the Hubbard model was first provided by Lieb and Wu [^], using 
Bethe Ansatz techniques. A particularly enlightening derivation is provided by Sutherland 



in Ref . fl^ . The resulting equations to determine the wavevectors of the N^. electrons are: 

2tt 1 ^ 

/cj = ^/i + ^E^(2sinA;,-2A^), j = 1, 2, 3, . . . , A^^. (3) 
An auxiliary set of M real numbers, the A^'s are determined through the set of equations 
E^(2sinfc,-2A,) = 27rJ, - ^^(A.-A^), a = 1, 2, 3, . . . , M, (4) 

j=l /3=1 

where 

9{p) = -2tan-\2pt/U), -tt < ^ < tt. (5) 

Here the Ij are consecutive integers (half-odd integers) if M is even (odd), and the Ja are 
consecutive integers (half-odd integers) if M' is odd (even). For the ground state we choose 
them to be (separately) clustered around zero (both negative and positive). In general, 
for finite size systems, their sum + J2aJa)) does not give zero, so that the total 

momentum of the ground state in some sectors is not zero. For the attractive Hubbard 
model M and M' are kept close to one another so that the maximum number of electrons 
with opposite spin can pair to take advantage of the attractive \U\. 

These equations can be iterated to convergence using a standard Newton-Raphson algo- 
rithm although some damping of solutions and judicious starting points are sometimes 
required. Once the kj have been obtained, the ground state energy is given by 

E = -2tJ2coskj. (6) 

One way to determine the energy of the first excited state is to compare the energy of two 
isolated systems with electrons in one and A^e + 2 in the other, with two isolated systems 
with Nf. + 1 in each. This yields the pair binding energy. 



efe(iV,) = 2E{Ne + 1) - E{Ne) - E{N, + 2). 



(7) 



Here we have used a single argument for the energy to indicate the total number of electrons 
— it is understood that M and M' are either equal (if A'^^; is even) or differ by one (if is 
odd). The single particle gap, A(Ai"e), is then given by 



In this model the gap value is positive for even and negative for odd, so that the overall 
curvature is in fact non-negative (see the second reference in 0). Thus, to calculate the gap 
value in this manner requires a distinction between even and odd numbers of electrons, a 
distinction which is unfortunately lost in the thermodynamic limit. Nonetheless, finite size 
effects can be monitored by increasing the system size until convergence is achieved. This 
becomes quite time-consuming, however, for weak coupling (where large system sizes are 
required for convergence), so that an alternative method is required, as discussed next. 

C. The Lieb-Wu Equations: Bulk Limit 

As oo, and A become distributed throughout the first Brillouin zone and the 

real axis, respectively, with density functions p{k) and cr(A), respectively ]12[. The above 
equations become the Lieb-Wu integral equations: 



A(A^e) = e(A^e)/2. 



(8) 





with cr(A) determined self-consistently from 



1 [Q \U\/At 



-p{k) dk = (t(A) + 



1 \U\/2t 



-ai\')dX', (10) 



71 J-Q {U/Atf + (A - sin k) 



TT J-B {U/2tf + (A - A') 



where Q and B are determined by the following sum rules: 




(11) 




(12) 
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with n = Ne/N and s = Sz/N. Finally the energy per lattice site is given by 

En cQ 

= -\U\{- - s) - 2t coskp{k)dk (13) 

i V ^ J — 

Note that we have already utilized Eq. (0) in these equations; hence they apply to an 
attractive Hubbard model with electron density n and magnetization s. The ground state 
is given by s = 0, so that Q = n, i.e. the entire Brillouin zone is occupied. These equations 
converge more easily than do their discrete counterparts given above, i.e. simple iteration 
is sufficient, and Newton- Raphson (and therefore the inversion of a large matrix) is not 
required. 

However, in principle, information concerning the gap is lost, as the distinction between 
even and odd numbers of electrons is no longer possible in the bulk limit. Nonetheless, 
Bahder and Woynarovich [jll| use the trick of applying an external magnetic field. Min- 
imization of the energy then determines the magnetization as a function of applied field. 
The magnetization is zero below a critical value of the field. The energy associated with 
this critical field gives the gap in the spin excitation spectrum. It turns out that half this 
energy corresponds exactly with the gap as defined by Eqs. (^ and (|]) . This is not a priori 
necessary, as the spin-gap is defined by flipping a spin with fixed electron number, whereas 
the previous gap was defined by changing the number of electrons in the system. 

To determine the gap in the bulk limit, we use a different procedure, which initially 
follows the method outlined in Ref. [ill] (^^^ their section IV where they carry out an 
analytical calculation in the low and high density limits). As already mentioned, the spin- 
gap can be defined as the energy required to flip one spin. Therefore we are interested in 
the energy difference between the s = state and the s = 1/N state where N is taken to 
approach oo. Since Q = tt for s = 0, we have from Eq. (|TT]) that Q = n — l/[A^p(7r)], using 
the even symmetry of p{k). Substitution into Eq. (^) yields one term on the left-hand side 
which is 0{1/N). Furthermore, p{k) can now be eliminated by substituting Eq. (^ into the 
remaining integral from — tt to vr [|ri|]: 
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^ r \U\/^t 2 \U\/At ^ IrB \U\/2t 

27r2 J-^ {U/Atf + (A - sin kf nN {U/Atf + ^ ' n J-b ^ ' {U/2tf + (A - A')^ ' 

(14) 

To first order in there is no change in the density function p{k)\ however, the same is 
not true for a(A). To see this we substitute s = 1/N into Eq. (p!2D , and allow B to adjust 
{B ^ Bo + SB) to the flipped spin (here the subscript o signifies the solution for s = 0). 
The important point is that the electron density n remains constant. Then 

where again we have utilized the even symmetry of cr(A). The two first terms within the 
braces do not cance/ because there are 0(1/ N) corrections to the density function cr(A). This 
is most clearly seen at half-filling, where 6B = and the sum rule has apparently changed 
by 0{1/N). Away from half-filling 6B is 0{1/N). Substitution of Eq. (|15|) into Eq. (|14D 
results in 



27r^J-^ {U/Aty + {X-smky n J-Bo ' ' {U/2ty + {X - X')^ 

2 \U\/At f \U\/2t \U\/2t ] 

TtN {U/Atf + A2 ^ TT ^ °^ I (U/2tY + (A - Bof {U/2tY + (A + B^Y J ^ ' 

where every term except the first on the right-hand-side contains (explicit or implicit) terms 

with 0{1/N). This equation is iterated to convergence (with high precision) with some 

choice for (eg. = 10^ or 10^). The equation for the energy, Eq. ( |l3|) becomes 

E nr.n 2t r ..2. \U\/^t 



'\U\--- dk cosH dX(x{X) 

2 TV J-TT J -Bo 



N ' '2 Ti J-^ J -Bo ' '(f//4t)2 + (A-sinA;)2 

, l^|-4t .^J 1^1/4 , 1^1/4 , 

N ^ °'\{U/At)^ + {Bo~smky {U/At)^ + {Bo + smk)^'' ^ ^ 



The gap energy is then identified as the coefficient of the 1/A^ term in the energy, 

E _ Eo 2A 
TV " A^ ^ aT' 



(18) 



The first term on the right-hand-side is obtained by solving the integral equations with 
A^ ^ oo, i.e. Eqs. (^0). Then the left-hand-side is obtained by solving Eqs. (|TB|-|I7p for 
A^ = 10^, 10^, etc., until the gap extracted from Eq. (|18|) has converged. 



D. The BCS Equations 



The derivation of the BCS equations is given in many places |14| and will not be presented 
here. They are: 

with 



E, = ^{e,-^^-\U\n/2f + ^lcs (20) 

and efc = —2t cos fc, where the sums are carried out over the first Brillouin zone (— vr < k < it). 
We use a finite lattice for simplicity, large enough for the gap to have converged to its bulk 
limit. An auxiliary equation to determine the electron number density is also required, 

1 ek- fi-\U\n/2 



For a given \U\ and n these equations are iterated to convergence to determine /x and Abcs- 
If the chemical potential lies within the band, then the gap to the first excited state, Ao, is 
given by Abcs- Otherwise, the gap is defined by the quasiparticle energy at the bottom of 
the band, i.e. 

Ao = Abcs for /i + |t/|n/2 > -2t 



(2t + |f/|n/2)2 + A|cs ioi ii+\U\n/2<-2t (22) 

Ao is to be directly compared to the gap discussed in the previous sections. Finally, the 
total energy is given by 

Note that the Hartree term is not required in the BCS gap equations (and is often omitted), 
but it provides an important contribution to the total energy. 

3. RESULTS 



A. The Ground State Energy 

In the dilute hmit the BCS solution of the gap equations is exact for all coupling strengths 
0. Of course the BCS approximation also reduces to the non-interacting result in the limit 
of weak coupling, for any electron density. Finally, BCS theory also becomes exact in the 
limit of strong coupling, as was appreciated for both the electron gas |[T^ and the lattice gas 



In Fig. 1 the ground state energy is shown as a function of electron density for various 
values of U. We have shown some results for repulsive U [0 along with the analytical results 
for U = 0{E/N = —^ sin ^) and U = oo {E/N = sinnvr). Note that agreement is best 
in weak coupling and in strong coupling, particularly for low electron density. Nonetheless, 
the BCS energy follows the exact result much more closely in the intermediate regime for 
this model than appears to be the case for the electron gas model 0. In fact a strong 
coupling expansion shows that the ground state energy has corrections of order [t/UY in 
both the exact and BCS solutions, whereas for the gas model the BCS solution contained a 
term linear in the inverse coupling constant which was not present in the exact solution . 

In Fig. 2 we show the ground state energy as a function of coupling strength, \U\ for 
half-filling {n = 1) and quarter-filling {n = 0.5). Here it is clear that the BCS theory is 
most inaccurate for intermediate coupling strength, i.e. \U\ ~ bandwidth. Nonetheless the 
maximum deviation is rather small, ~ 4%, which occurs at half-filling near \U\ = At. Thus, 
on the basis of ground state energy calculations, BCS theory appears to be a very accurate 
theory, even in one dimension. However, in the next subsection we show that the BCS result 
for the energy gap is far less accurate. 

B. The Energy Gap 

The calculation of the energy gap is a much more sensitive test of the accuracy of BCS 
theory. In Fig. 3 we show the g function of electron density for various values of the 

interaction strength, \U\. As already stated, the BCS gap is exact in the low density limit, 
regardless of the coupling strength. However, as the electron density increases, the true 



gap decreases for all coupling strengths, while the BCS gap (Ao) from Eq. ( p2D increases 
for all coupling strengths. The true gap decreases monotonically as a function of electron 
density, whereas the BCS gap is in general non-monotonic, exhibiting a maximum at some 
intermediate electron density, which is a function of coupling strength. For sufficiently 
strong coupling, the BCS gap is monotonically increasing, so that the maximum occurs at 
half-filling. The most serious errors occur near half-filling, particularly in the weak coupling 
limit, where BCS theory vastly overestimates the value of the gap. 

At half-filling certain analytical results can be obtained. For the true gap we obtain 

where Ji{uj) is the Bessel function of the first kind of order 1. This result can be written in 
a more useful form ||18||, 



where Ki{x) is the first order modified Bessel function. In weak coupling this reduces to an 
exponential contribution: 

-27rt 



Aweak = -Vl?7|texp|^^[>. (26) 



In strong coupling the first two terms of Eq. ( p4D dominate. 

The BCS integral can also be performed analytically: we obtain, at half-filling, 

2vrt 1 / 1 \ 

where 6 = ^ and K{k) = Jq^'^ . ^ =d9 is the complete elliptic integral of the first kind. 
Note that Ao = Abcs at half-filling. In weak coupling we obtain 

Ao(weak) = 8texp(^-^^ (28) 
while in strong coupling we have 

A _ M _ ^ 

^o{strong) — 2 • V^/ 
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Thus, the weak couphng BCS gap at half-filhng has the famihar activated form, albeit with 
incorrect prefactors. In strong couphng the first order correction (of order unity) is absent 
in the BCS solution. 

Finally, at zero filling, the exact solution is known analytically |T7| : 



A = -2t+^(M)V(2tf, (30) 

and the BCS result is identical. Note that in weak coupling the dependence on \U\ is 
quadratic. 

In Fig. 4 we show the gap as a function of coupling strength for zero-, quarter-, and 
half-filling. As already stated, for zero filling BCS is exact, while for any non-zero filling, the 
deviations are as shown. BCS theory overestimates the gap by a considerable degree over the 
entire coupling range shown. Note that the variation of the exact gap with electron density 
is not significant beyond \U\/t ^ 2, and, as remarked earlier, monotonically decreases with 
increasing electron density. The dependence of the BCS gap on the electron density is not 
monotonic, and beyond \U\/t ^ 3 the ordering is incorrect. In Fig. 4b we show the weak 
coupling regime in greater detail. In this regime the ordering with electron density is correct, 
but BCS theory greatly overestimates the gap value, as shown. 

C. Finite Size Effects 

The attractive Hubbard model can be solved in the thermodynamic limit, and therefore, 
for this study, finite size effects do not concern us. However, most models require some sort 
of numerical solution, often on small system sizes, and are therefore subject to finite size 
effects. It is therefore of interest to examine how these effects influence the results in both 
the exact and BCS solution in the present model. To some extent the following conclusions 
can be applied to other models where an exact solution in the bulk limit is not possible. 

To obtain the exact gap it is necessary to solve Eqs.(0-|^) for some fixed lattice size, 
N. The BCS gap is obtained from Eqs.(^^-^^, again for a fixed lattice size, A^. However, 



these equations utilize a grand canonical ensemble, so that particle number is fixed only in 
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an average sense, through the chemical potential. Thus, all electron densities are possible 
(except for certain densities where discontinuities may arise); this raises the issue of how 
applicable the use of the grand canonical ensemble is for small clusters. We partially address 
this issue here by using the canonical ensemble for clusters with one, two, three and four 



electrons. This has already been done for a tetrahedral cluster [113, where, for symmetry 
reasons, the BCS solution is exact for all pair-fillings on a four-site cluster. We use the BCS 
wavefunctions 



^4> = E C^(/^i,A;2)4^cLfc,i4T^UilO> (31) 



fe2>fel 



for two and four electrons, respectively. On a chain of four sites these would correspond 
to quarter- and half-filling. To accomodate odd numbers of electrons we simply add an- 
other electron, whose momentum then determines the total momentum of the ground state 
wavefunction. For example, for three electrons, 

I*3>,= E^?^4tcVJtI0> (32) 

k^q 

and the third electron has the important role of blocking the pair from adopting the (g, —q) 
pair state. Minimization with respect to the relevant parameters yields equations for the 
or C{ki, k2) from which the minimum energy can be obtained. The gap energy can then be 
obtained following the method outlined in Section 2.B. In contrast to the grand canonical 
method for finite systems, this method generates gap energies for values of the electron 
density that are commensurate with the lattice size. 

In Fig. 5 results for the gap are illustrated as a function of density for \U\ = — 2t, using 
all the various methods discussed thus far. The open symbols represent exact solutions 
for various lattice sizes as indicated in the figure caption. The gap generally decreases as 
the lattice approaches the bulk limit. The same trend is observed in the BCS solutions, 
indicated by the various curves. Finally, for = 4, the filled circle indicates the result for 
quarter filling from the canonical BCS equations. 
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We have found quite generally that the use of the canonical ensemble always improves the 
ground state energy, compared to that obtained from the grand canonical ensemble. This 
is particularly true for electron densities that correspond to an odd number of electrons, 
and for low electron densities. In fact the conventional BCS theory is variational with 
respect to H — /iN, and not the Hamiltonian alone, so that for small system sizes this 
theory gives lower energies than the true ground state energies at low electron densities. 
This does not violate the variational principle since the BCS wavefunction used spans all 
particle numbers, and BCS theory never gives a lower energy than the true ground state 
energy for any particle number. In any event, this distasteful feature is remedied by the use 
of the canonical ensemble. Of course, as the bulk limit is approached these discrepancies 
disappear. 

However, the use of the canonical ensemble does not universally improve the BCS gap 
energy for small systems. The case shown in Fig. 5 shows some improvement. For weaker 
coupling the improvement is insignificant, while for stronger coupling the result from the 
canonical solution is poorer than the grand canonical one. 

It should be clear that finite size effects are significant only in the intermediate to weak 
coupling regime (the two site system reproduces the correct bulk limit in the strong coupling 
limit). For the intermediate coupling case shown, the finite size results could mislead one to 
believe that the gap energy increases as electron density increases from zero. This behaviour 
is displayed by the BCS results as well, and remains a feature of the BCS solution in the 
bulk limit (solid line). However, the true energy gap is a monotonically decreasing function 
of electron density. Furthermore, the gap always decreases as the bulk limit is approached, 
so that the relative error obtained from a finite system study will always underestimate the 
true error in the bulk limit. The error in the BCS gap at half-filling, for example, is 20% for 
a four site chain, while, in the bulk limit it is 300%. This underestimate worsens for weaker 
coupling. 

Finally, we should note that the finite size effects on the ground state energy are not 
nearly as severe as just discussed for the energy gap. For example, the ground state energy 
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at half- filling for N — 8 differs from the bulk limit by about 1% while the gap value is almost 
a factor of three too large. 

4. SUMMARY 

We have carried out an evaluation of BCS theory for the attractive Hubbard model in 
one dimension. This is certainly the simplest Hamiltonian with many-body interactions for 
which an exact solution exists, so that an evaluation in the bulk limit is possible. We hope 
that the results of this paper can be used for more complicated Hamiltonians where exact 
solutions are not possible. 

The ground state energy is very accurately reproduced by the BCS wavefunction, par- 
ticularly at low electron densities, and for all coupling strengths. In fact the ground state 
energy is least accurate at intermediate coupling strengths, at half-filling, and here the error 
is only a few per cent. This result is encouraging, and ought to improve with increasing 
dimensionality. 

The energy gap is not so accurately reproduced within BCS theory. There is exact 
agreement in the dilute limit; however the exact result monotonically decreases with electron 
density while the BCS result increases, at least initially. There is also agreement in the strong 
coupling limit, and BCS theory displays the same activated form in the weak coupling 
limit, but with an incorrect prefactor. At intermediate coupling strength the BCS gap 
overestimates the true gap typically by a factor of two or more. The relative error decreases 
(increases) with stronger (weaker) coupling. Over most of the coupling strength regime the 
maximum error occurs at half-filling. 

Finally, wc have taken advantage of this exactly solvable case to examine finite size 
effects. Insofar as these effects occur in other more complicated models, they can be used 
as cautionary guidelines in future finite system studies. 

The relevance to higher dimensionality remains somewhat of an open question. While 
the accuracy of the BCS ground state energy is most encouraging, it is also clear from this 
study that this agreement is not a good indicator of the accuracy of the BCS energy gap. 
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On the other hand the poor accuracy of the BCS gap in one dimension may well be rectified 
by proceeding to higher dimension. It is hoped that future investigations of this issue can 
utilize this study cis db benchmark. 
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FIGURES 



Fig. 1 Marsiglio 




FIG. 1. The ground state energy as a function of electron density for various values of the 
coupling strength, both repulsive {U > 0) and attractive {U < 0). BCS results are also shown for 
comparison. 
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Fig. 2 Marsiglio 
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Fig. 3a Marsiglio 
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FIG. 3. The energy gap, A/i vs. electron density for (a) \U\/t = I, (b) \U\/t = 2, and (c) 
\U\/t = 10. The BCS result (i.e. Ao/t) is shown for comparison. 
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Fig. 3c Marsiglio 
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Fig. 4a Marsiglio 




FIG. 4. (a) The energy gap, A/f vs. coupling strength, \U\/t, for various electron densities. In 
(b) we provide an expanded version of the weak coupling regime. 
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Fig. 4b Marsiglio 
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Fig. 5 Marsiglio 




FIG. 5. The energy gap, A/f vs. electron density for \U\/t = 2, for various chain sizes. The 
open symbols are from exact results, with lattice size as indicated. The curves are BCS results 
(using the grand canonical ensemble) with lattice sizes as indicated. We have also indicated the 
result for quarter-filling for N = 4 (filled circle), showing an improvement over the grand canonical 
result. 
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